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Research  article 

Phase  Field  Theory  and  Analysis  of  Pressure-Shear  Induced 
Amorphization  and  Failure  in  Boron  Carbide  Ceramic 


John  D,  Clayton 

Impact  Physics  Branch,  US  Army  Research  Laboratory,  Aberdeen  MD  21005-5066,  USA;  Email: 
john.d.elaytonl.eiv@mail.mil;  Tel:  +1-410-278-6146;  Fax:  +1-410-278-2460. 

Abstract:  A  nonlinear  eontinuum  phase  field  theory  is  developed  to  deseribe  amorphization  of 
erystalline  elastie  solids  under  shear  and/or  pressure  loading.  An  order  parameter  deseribes  the  loeal 
degree  of  erystallinity.  Elastie  eoeflheients  ean  depend  on  the  order  parameter,  inelastie  volume 
ehange  may  aeeompany  the  transition  from  erystal  to  amorphous  phase,  and  transitional  regions 
parallel  to  bands  of  amorphous  material  are  penalized  by  interfaeial  surfaee  energy.  Analytieal  and 
simple  numerieal  solutions  are  obtained  for  an  idealized  isotropie  version  of  the  general  theory,  for 
an  element  of  material  subjeeted  to  eompressive  and/or  shear  loading.  Solutions  eompare  favorably 
with  experimental  evidenee  and  atomie  simulations  of  amorphization  in  boron  earbide, 
demonstrating  the  tendeney  for  struetural  eollapse  and  strength  loss  with  inereasing  shear 
deformation  and  superposed  pressure. 

Keywords:  eeramies;  phase  transformations;  amorphization;  boron  earbide;  phase  field;  elastieity 


1.  Introduction 

Many  crystalline  ceramics  and  minerals  undergo  structural  changes  (e.g.,  phase  transformations, 
twinning,  or  fracture)  when  subjected  to  stresses  of  high  magnitude.  One  such  structural  transition  is 
amorphization,  i.e.,  transformation  from  a  crystalline  structure  to  an  amorphous/glassy  phase  lacking 
long-range  order.  Under  mechanical  loading  by  large  compressive  and/or  shear  stress,  such 
amorphization  may  occur,  for  example,  in  quartz  (a-SiOi;  rhombohedral  crystal  structure)  [1], 
berlinite  (a-AlP04;  rhombohedral  structure)  [2],  some  forms  of  garnet  (Y3AI5O12;  cubic  structure) 
[3],  and  boron  carbide  (nominally  B4C,  with  varying  carbon  content  possible;  rhombohedral)  [4]. 
Noteworthy  among  all  of  these  materials  is  a  loss  of  static  or  dynamic  shear  stiffness  with  increasing 
pressure.  According  to  a  theoretical  criterion  attributed  first  to  Bom  [5]  and  advanced  in  subsequent 
literature  [6-8],  attainment  of  null  stiffness  in  one  or  more  preferred  directions,  which  can  be 
associated  with  loss  of  positive-definiteness  of  a  certain  tangent  elastic  stiffness  tensor  depending  on 
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loading  protocol,  may  signal  onset  of  a  localized  structural  transformation  such  as  amorphization. 

The  present  paper  focuses  on  boron  carbide,  a  strong  ceramic  of  high  stiffness,  high  hardness, 
low  mass  density,  and  low  ductility.  Different  structural  variants  of  boron  carbide  exist,  called 
polytypes.  Atomic  chains  consisting  of  C-C-C  or  C-B-C  are  aligned  parallel  to  the  [0001]  direction 
(hexagonal  Miller  indices)  and  thread  basal  layers  of  icosahedra  at  the  rhombohderal  vertices  of  the 
unit  cell.  Polytypes  may  differ  in  stoichiometry,  ground  state  energy,  elastic  stiffness,  and 
amorphization  tendency  [9-11].  One  model  for  amorphization  based  on  results  of  Density  Functional 
Theory  (DFT)  calculations  suggests  that  structural  collapse  and  localization  follow  cross-linking  of 
the  C-B-C  chain  with  icosahedral  atoms  [4];  other  calculations  suggest  segregation  is  most  favorable 
for  polytypes  with  C-C-C  chains  [9].  More  recent  DFT  calculations  [12]  demonstrate  amorphization 
in  both  classes  of  polytype  under  uniaxial  loading.  In  physical  experiments,  amorphization  of  boron 
carbide  has  been  observed  in  Diamond  Anvil  Cell  (DAC)  compression-decompression  [4], 
indentation  [13],  and  ballistic  impact  [14].  Regarding  the  latter,  performance  of  boron  carbide 
ceramic  in  protection  systems  such  as  armor  for  vehicles  and  personnel  is  thought  to  be  severely 
impeded  by  its  tendency  to  localize,  with  cleavage  fracture  accompanying  or  closely  following 
amorphization  [14].  In  recovered  samples,  regions  of  glassy  phase  are  often  reported  to  be  in  the 
form  of  planar  bands  of  small  thickness,  on  the  order  of  several  nanometers,  and  may  be 
preferentially  located  parallel  to  certain  crystallographic  planes  [4,  14]. 

In  order  to  facilitate  design  of  structures  and  composite  material  systems  for  defense  and 
industrial  applications,  continuum  scale  models  are  needed  for  predicting  localization  and  failure  in 
boron  carbide  ceramics.  Previous  models  for  the  transition  from  crystalline  to  glassy  phase  include 
one  based  on  adiabatic  shear  localization  [15]  and  another  in  which  instability  is  triggered  by  loss  of 
nonlinear  elastic  tangent  stiffness  [16,  17]  similar  to  Born’s  criterion.  In  contrast,  the  present  paper 
invokes  nonlinear  continuum  phase  field  theory,  in  what  appears  to  be  the  first  known  application  of 
phase  field  modeling  towards  stress-induced  amorphization. 

In  the  phase  field  approach,  the  thermodynamic  state  of  the  material  is  described,  locally  and  in 
part,  by  one  or  more  order  parameters.  Elastic  strain  energy  density  generally  depends  on  order 
parameter(s),  as  does  surface  energy,  the  latter  reflected  by  dependence  of  total  energy  on  spatial 
gradients  of  the  order  parameter(s).  Phase  field  models  have  been  applied  to  describe  numerous 
kinds  of  structural  changes  in  crystalline  materials,  including  dislocation  glide  [18],  twinning  [19-21], 
and  fracture  [22-25].  Advantages  of  phase  field  methods  over  other  continuum  models  for  crystal 
plasticity,  twinning,  and  fracture  in  metals  [26]  and  ceramics  [27-29]  include  the  following:  (i) 
relatively  few  material  parameters  are  needed,  (ii)  evolution  laws  for  inelastic  deformation 
mechanisms  follow  naturally  from  energy  minimization  and  need  not  be  prescribed  ad  hoc,  and  (iii) 
an  intrinsic  length  scale  is  introduced  via  the  surface  energy,  enabling  regularization  of  localization 
zones  and  mesh  independent  numerical  results.  Regarding  the  latter  advantage,  the  need  for  explicit 
tracking  of  sharp  interfaces  is  avoided  in  numerical  schemes,  and  standard  finite  element 
discretizations  may  be  used  for  solution  of  boundary  value  problems. 

In  the  present  paper,  phase  field  theory  is  applied  towards  amorphization  in  a  way  similar  to 
prior  phase  field  models  for  fracture  [22-25],  with  a  single  order  parameter  increasing  in  value  from 
zero  to  unity  during  the  transition  from  crystalline  to  amorphous  phase,  and  with  shear  stiffness 
degrading  with  increasing  order  parameter.  It  is  clarified  that  whereas  amorphization  can  be 
classified  as  a  true  phase  transformation,  twinning  and  fracture  are  not  usually  categorized  as  such  in 
a  strict  sense.  Regardless,  phase  field  models  tend  to  treat  twinning  and  fracture  analogously  to  phase 
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changes,  wherein  evolution  of  order  parameter(s)  quantify  the  transition  in  state  from  a  perfect  parent 
erystal  “phase”  to  a  fully  twinned  or  fully  fraetured/damaged  “phase”.  Struetural  collapse  may  be 
aeeompanied  by  an  inelastie  decrease  in  volume  (i.e.,  increase  in  mass  density)  [4],  which  promotes 
amorphization  with  inereasing  pressure.  Thus,  under  pressure,  shear  stiffness  and  strength  tend  to 
deerease,  in  agreement  with  DFT  calculations  [10-12].  In  the  present  paper,  a  general  thermodynamic 
theory  is  developed  first,  aceounting  for  anisotropic  elastic  strain  energy  and  anisotropie  surface 
energy.  Then  an  idealized  nonlinear  isotropic  theory  is  developed  that  facilitates  solution  of  boundary 
value  problems.  Specifieally,  problems  involving  volumetric  compression  and  simple  shear  and  are 
studied  analytieally  to  demonstrate  the  predietive  ability  of  the  theory.  Stability  is  assessed, 
extending  previous  analyses  of  eontinuum  damage  models  [30].  Effects  of  properties  particular  and 
peculiar  to  boron  carbide  entering  the  model  are  explored,  providing  new  insight  into  the  onset  of 
loealization  and  failure.  Notation  of  nonlinear  eontinuum  meehanics  [31]  is  used,  here  primarily  in 
direct  form  for  vectors  and  tensors  which  are  written  in  bold  font. 

2,  Materials  and  Method 

2.1.  Boron  Carbide 

Single  crystals  of  boron  carbide  (B4C)  belong  to  space  group  R3m,  centrosymmetrie  point 
group  3m,  and  Laue  group  RI,  the  latter  indieative  of  six  independent  seeond-order  elastic  constants 
and  fourteen  independent  third-order  elastic  constants  [16,  31].  The  theoretical  mass  density  of 
stoiehiometric  B4C  is  2.52  g-cm'  .  Polyerystalline  boron  carbide  has  a  typical  grain  size  on  the  order 
of  10  |um,  an  elastic  modulus  on  the  order  of  470  GPa,  a  hardness  on  the  order  of  30  GPa,  and  a 
Hugoniot  Elastic  Eimit  (HEE,  or  yield  stress  under  shoek  compression)  on  the  order  of  15-20  GPa 
[17].  Ductility  is  low,  implying  that  although  full  dislocations,  partial  dislocations,  and  stacking 
faults  have  been  observed  [32],  their  low-temperature  mobility  is  inhibited  [16].  Twins  are  thought  to 
arise  predominantly  during  processing  [33]  rather  than  to  be  induced  by  stress  or  strain.  Primary 
inelastic  deformation  meehanisms  under  impact  loading  are  amorphization  and  fracture  (primarily 
cleavage  fraeture  on  intrinsically  weakest  planes),  which  may  occur  sequentially  or  simultaneously. 


Table  1.  Properties  of  boron  carbide. 


Property 

Description 

Value  [units] 

Bo 

bulk  modulus  (crystal) 

248  GPa 

Po 

shear  modulus  (crystal) 

197  GPa 

/ 

localization  width  (glass) 

2  nm 

r 

surface  energy  (glass) 

347  mJ-m'^ 

M 

density  ratio  (crystal/glass) 

0.96 

Eisted  for  future  reference  in  Table  1  are  physieal  properties  entering  the  idealized  isotropie 
phase  field  model  of  boron  carbide  to  be  applied  later  in  analysis  of  pressure-shear  deformation  and 
eorresponding  instability.  Ambient  bulk  modulus  Bq  and  shear  modulus  po  for  the  fully  erystalline 
phase  are  obtained,  respectively,  from  a  fit  of  the  nonlinear  elastic  potential  function  to  quantum 
mechanical  equation-of-state  data  for  the  polar  C-B-C  polytype  [11]  and  a  Voigt  average  of  the 
second-order  anisotropic  constants  [29,  31].  Regularization  parameter  /  indicating  the  width  of 
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localized  structural  transition  zones  is  representative  of  the  thiekness  of  planar  glassy  phases  of 
boron  earbide  observed  in  recovered  experimental  samples  [4,  14],  The  surfaee  energy  of  amorphous 
zones  is  computed  as  the  difference  between  ground  state  energies  (measured  per  appropriate  unit 
reference  volume)  of  segregated  amorphous  and  crystal  phases,  multiplied  by  /,  where  energy  of  the 
polar  C-B-C  polytype  is  used  for  the  crystal  reference  [9].  The  ratio  of  mass  density  of  the  erystalline 
phase  to  that  of  the  amorphous  phase  is  denoted  by  S;  this  value  depends  on  the  composition  and 
structure  of  the  glass,  here  taken  as  0.96,  corresponding  to  a  4%  volume  reduction  upon  structure 
collapse  eommensurate  with  amorphization  [4], 

2.2.  Phase  Field  Theory:  General 

General  nonlinear,  anisotropic,  three-dimensional  (3D)  theory  is  outlined  in  §2.2.  Model 
kinematies,  energy  functions,  and  governing  equations  are  described  in  turn. 

2.2.1.  Kinematies 

Let  x=x  {X,t)  denote  spatial  position  at  time  t  of  a  material  particle  whose  reference  position  is 
denoted  by  X.  The  deformation  gradient  is 

F=Vx  =  dxldX  (1) 

where  V  is  the  material  gradient  operator.  Let  [0,1]  denote  the  scalar  order  parameter,  where 
=  0  V2f  e  crystal  phase  at  time  t 

^(2f ,  t)  <  e  (0, 1)  V2f  e  interface  zone  at  time  t  (2) 

=  1  \/X  e  glassy  phase  at  time  t 

The  deformation  gradient  is  split  multiplieatively  into  a  product  of  two  terms  as 

F  =  F^F\  F\X,t)  =  F\i,{X,t)-\  (3) 

where  is  the  elastic  deformation  conjugate  to  the  applied  stress  and  is  the  inelastic  deformation 
that  depends,  as  indieated,  on  the  loeal  value  of  the  order  parameter  associated  with  structural 
transformation,  specifically  amorphization.  Unlike  F,  whieh  always  obeys  compatibility  (null  curl) 
conditions  V  x  F  =  0,  neither  F^  nor  F^  is  individually  always  integrable  to  a  vector  field,  i.e.,  these 
two  deformations  are  generally  anholonomic  [34].  Note  that  this  theory  only  accounts  explicitly  for 
one  type  of  structural  change.  Not  included  in  (3)  are  additional  terms  needed  for  representation  of 
effects  of  other  kinds  of  defects,  for  example  disloeation  slip  [26,  31],  twinning  [27],  fraeture  [35], 
diselination  rotation  [36],  and  dilatation  assoeiated  with  stacking  faults  and  dislocation  or 
disclination  cores  [37,  38]  or  point  defects  [39,  40].  The  symmetric  elastic  deformation  tensor  C  and 
the  measure  of  elastic  volume  ehange  J  are  defined  as  (no  “E”  superseript) 

C  =  (F^)V^  J  =  Vdet  C  =  det  =  det  F/det  (4) 

2.2.2.  Energy  Eunctional 

Similar  to  previous  nonlinear  phase  field  theories  for  twinning  [19-21]  and  fracture  [25],  a  total 
energy  functional  for  a  body  of  reference  volume  V  is  posited  in  general  form  as 
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^(x,  =  I  [WiC,  +  m  +  K :  (X)  V^W  (5) 

where  W  is  elastic  strain  energy  density  per  unit  reference  volume,  /  accounts  for  the  difference  in 
ground  state  energies  among  crystalline,  amorphous,  and  intermediate  phases,  and  k  is  a  symmetric 
second-order  tensor  that  accounts  for  surface  energy  associated  with  gradients  of  the  order  parameter. 
For  a  fully  realistic  nonlinear  anisotropic  response,  W  should  reflect  the  rhombohedral  symmetry  of 
the  crystalline  phase,  the  importance  of  higher-order  elastic  constants  [16],  and  the  transition  to 
elastic  isotropy  upon  complete  amorphization  [17].  In  order  to  restrict  amorphous  bands  to  specific 
crystallographic  planes,  the  following  form  of  K  may  be  used  [25]; 

K  =  Ko[1-i-[3(1-»i  0/m)]  (6) 

Here  kq  is  a  scalar  and  [3  >  0  is  a  penalty  factor  that,  when  greater  than  zero,  increases  surface 
energies  on  planes  not  normal  to  referential  unit  vector  mi.  If  such  a  representation  is  insufficient  (i.e., 
if  numerous  different  planes  with  various  surface  energies  and  orientations  are  prone  to  localization), 
the  theory  may  require  extension  to  distinct  order  parameters  for  each  plane. 

2.2.3.  Governing  Equations 

A  variational  framework  is  considered  henceforth,  amenable  to  solution  via  energy  minimization 
[19-21].  Time  t  becomes  a  parameter  and  does  not  explicitly  enter  the  governing  equations,  which 
are  now  applicable  only  to  quasi-static  problems.  The  following  variational  principle  is  set  forth; 

5'F  =  |jG5jc  +  r5^]d5  (7) 

with  t  the  mechanical  traction  vector,  r  a  conjugate  force  to  the  order  parameter,  and  S  the  referential 
surface  enclosing  V.  Application  of  standard  mathematical  procedures  involving  integration  by  parts 
and  the  divergence  theorem  [19,  31]  results  in  the  following  local  equilibrium  equations  (i.e., 
Euler-Lagrange  equations)  in  V  and  natural  boundary  conditions  on  S: 


VP  =  0, 

P  =  {dW/dF)[^ 

(8) 

q  +  d//d^  =  2KV^^, 

?=(e»f/05)|. 

(9) 

t  =  Pn, 

r  =  2k;  V^0M 

(10) 

The  first  Piola-KirchhofT  stress  tensor  is  P,  and  the  unit  outward  normal  vector  to  S'  is  m. 

2.3.  Phase  Field  Theory:  Idealized 

The  nonlinear  theory  of  §2.2  is  now  specialized  to  isotropic  behavior.  This  model,  which  is  far 
more  amenable  to  analytical  and  simple  numerical  solutions  than  that  of  §2.2,  can  be  applied  to  gain 
insight  into  particular  problems  (e.g.,  ID)  wherein  anisotropy  is  of  secondary  importance.  Model 
kinematics,  energy  functions,  and  governing  equations  are  again  described  in  turn. 
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2.3.1.  Kinematics 

Equations  (l)-(4)  still  apply.  Inelastic  deformation  corresponding  to  amorphization  is  specified 
as  the  isotropic  volume  change  (i.e.,  spherical  eigenstrain) 

F'©  =  <K0)  =  1,  <K1)  =  3  (11) 

with  (j)  a  sufficiently  smooth,  positive  scalar  function  of  the  order  parameter  and  with  5  introduced  in 
§2.1.  A  similar  construction  has  been  used  elsewhere  to  model  residual  volume  changes  from  point 
defects  [40]  and  pore  collapse  [41].  Equation  (4)  then  reduces  to 

C(F,0  =  [^(0]-^'^F V  ,  =  detF/(^(^)  (12) 

2.3.2.  Energy  Eunctional 

Equations  (5)  and  (6)  still  apply.  A  compressible  neo-Hookean  nonlinear  elastic  strain  energy 
potential  [20,  21,  25]  is  prescribed; 

W{C,i,)  =  \\i{ixC -3 -2\nJ)  +  \X{\njf  (13) 

Elastic  coefficients  may  depend  on  the  order  parameter 

E  =  |LioX(^),  =  ^o=^o-|Eo’  X(0)  =  i(0)  =  l  (14) 

Here,  y  and  i  are  sufficiently  smooth  scalar  functions  of  the  order  parameter  that  interpolate  between 
elastic  constants  of  crystal  and  glass  phases.  Isotropic  bulk  and  surface  energies  of  the  amorphous 
phase  are  represented  by 

/(0  =  r^7/,  K;v^0v^  =  r/|v^f  (K„=r/,p  =  o)  (15) 

where  surface  energy  per  unit  reference  area  E  and  regularization  length  /  have  been  introduced  in 
§2.1.  Total  energy  (5)  then  becomes 

'E(x,^)=f  [ffi(Vx,^)  +  r7//  +  r/|V^f  ]dE  (16) 

with  strain  energy  (13)  expressed  explicitly  in  terms  of  ^  and  F=Vx  as 

ffi(F,^)  =  iqoX(0{[^(07'^F:F-3-21n[detF/(^(0]}  +  iXoi(^){ln[detF/^(0]f 

(17) 

As  a  point  of  clarification,  the  term  “interpolation  function”  is  used  herein  to  refer  to  any  of  functions 
(j),  X,  and  i  that  depend  continuously  on  values  of  order  parameter  These  functions  are  used  to 
estimate  values  of  specific  volume  and  elastic  coefficients  in  material  wherein  transformation  is 
incomplete  (i.e.,  wherein  0<^<  1),  given,  a  priori,  known  values  of  such  quantities  at  the  endpoints 
^  =  0  and  ^  =  1  corresponding  to  pure  crystalline  and  amorphous  phases,  respectively. 

2.3.3.  Governing  Equations 

Equations  (7)-(10)  still  apply.  Using  (13),  stress  in  the  second  of  (8)  and  elastic  driving  force  in 
the  second  of  (9)  become,  respectively. 
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P{F,  0  F-F-^}  +  ln[det  F/mW^"  ( 1 S) 

q{F,  0  ^  '(0  {^0  X(0/^(0  -  i  F  .F-  Xo [i(0/^(0]  ln[det  F/^(0] } 

+  X  li  ^ ^  - 1  -  ln[det  F/^(0] }  +  i  i  {ln[det 

(19) 

where  (j)'  =  d(t)/d^,  x'  =  dx/d^,  and  i'  =  di/d^.  Cauehy  stress  tensor  a  =  (l/detF)PF^  and  Cauehy 
pressure  p  =  -tra/S  are,  respeetively, 

a(F,  0  =  [^ox(0/det  F]  { FF"  - 1}  +  [X^iCO/det  F]  ln[det  F/^(0]1  (20) 

p{F,  0  =  [^ox(0/det  F]  (1  -  i  [^(0]-^'^  F  ;  F}  -  [Xoi(^)/det  F]  ln[det  F/^(0]  (2 1 ) 


3,  Results 

Analyzed  are  solutions  of  the  phase  field  theory  of  §2.3  for  spatially  homogeneous  hydrostatie 
eompression,  followed  by  solutions  for  eombined  spherieal  eompression  and  simple  shear. 

3.1.  Hydrostatic  Compression 

First  eonsider  uniform  spherieal  deformation  of  an  element  of  initially  fully  erystalline  material 


from  referenee/initial  volume  V  to  eurrent  volume  v: 

x  =  A''^X;  F  =  A''^1;  ^  =  detF=v/F>0  (V^  =  0)  (22) 

Under  sueh  spherieal  deformation,  (18)-(21)  reduee  to 

F(A,^) = nmim  (23) 

q{A,^)  =f(^)K  xm(^)-fioX(0[c^(^)r''^'''  -X„[i(^)/(^(^)]ln[z(/(^(^)]} 

+  X  ’(^)fio  {f  -  f  -  HA/m] }  +  i  I  ’(^)Xo  {ln[z(/(^(^)] 

cy{A,^)  =  M^)lA-\  { [^(0]-^'^  - 1}1  +  [Xoi(0/^]  ln[^/^(0]l  (25) 

p{A,^)  =  [FoX(0/^]  (1  -  [Almr  }  -  [^oi(^)/^]  HAlm^  (26) 

In  Cartesian  eoordinates  x  =  {x,y,z)  and  X  =  {X,Y,Z),  stress  equilibrium  in  the  first  of  (8)  with  (22) 
and  (23)  reduees  to 

dP^ldX  =  dP^.jdY  =  dP^JdZ  =  t)^P  =  P^^=  P^^  =  =  constant  (27) 

VPiA,  ^)  =  idP/dA)VA  +  (0F/0^)V^  =  0  ^  ^  =  constant  (28) 


Order  parameter  equilibrium  in  the  first  of  (9)  with  equations  (15),  (24),  and  (28)  becomes  the 
algebraic  equation 
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^ = -  i  (//r)^  ’[^0  x/^  -  -  ^0  ( v^)  nm)] 

-K//r)xVo[f^^'''^'''-f-ln(^/^)]-i(//r)i'Xo[ln(^/^)f 

Solution  of  this  equation  requires  speeifieation  of  interpolation  funetions  (j)(^),  x(^),  and  i(^).  Similar 
to  the  polynomial  form  used  elsewhere  in  a  phase  field  model  for  twinning  [19-21],  let 

^  =  S  +  (1  -  S)[3(l  -  -  2(1  -  f  =  -6(1  -  S)[(l  -  0  -  (1  -  ]  (30) 

Similar  to  phase  field  models  for  fraeture  invoked  elsewhere  [22-25],  the  shear  modulus  degrades  as 

X  =  (l-^)^  X'  =  -2(1-^)  (31) 

This  is  eonsistent  with  the  reported  loss  of  shear  strength  in  loealized  amorphous  bands  in  boron 
earbide  [4,  14,  16].  Transformation  ean  also  be  interpreted  to  eorrelate  with  mode  II  fraeture  in  this 
model.  Now  let  the  ambient  bulk  modulus  5o  remain  eonstant  upon  amorphization,  whieh  requires 
from  (31)  that 

1  =  - 1(1  -  n.  A,  =  1(1  +  V.  -  (1  -  5)=(1  -  2v.)]/v. 

i'  =  |(l-5)(l-2v.)/v„ 

with  Vo  Poisson’s  ratio  for  the  erystalline  phase.  Prior  DFT  ealeulations  [42]  have  suggested  a  w4% 
deerease  in  bulk  modulus  may  oeeur  upon  strueture  transformation;  this  ehange  is  eonsidered  small 
enough,  and  of  the  same  order  of  aeeuraey  as  the  DFT  solutions,  to  be  justifiably  ignored,  thereby 
simplifying  the  theoretieal  framework  and  analytieal  solutions. 


Figure  1.  Pressure-volume  response  (a)  and  order  parameter  evolution  (b)  under 
uniform  spherical  compression. 

Solution  data  are  obtained  as  follows  for  hydrostatie  eompressive  loading.  Volume  is  redueed 
inerementally  by  deereasing  A  from  an  initial  eondition  of  unity.  For  eaeh  volume  deerement,  (29)  is 
solved  via  simple  numerieal  iteration  to  obtain  order  parameter  Onee  ^  is  known,  funetions 
(30)-(32)  are  evaluated  and  then  substituted  into  analytieal  expression  (26)  to  determine  pressure  p. 
Compared  in  Figure  1(a)  are  the  pressure  predieted  from  the  present  phase  field  model  and  results  of 
quantum  (DFT)  ealeulations  [11]  for  the  polar  C-B-C  polytype,  whieh  were  fit  to  a  Bireh-Murnaghan 
equation-of-state  [43].  Agreement  between  the  two  sets  of  results  is  deemed  exeellent  for  ^  >  0.9,  i.e., 
for  volumetrie  eompression  less  than  around  10%.  As  shown  in  Figure  1(b),  ^  remains  very  small  for 
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compressions  up  to  15%;  in  fact,  the  phase  field  model  predicts  0.1  for  A  >  0.64.  Such 
predictions  agree  with  DFT  simulations  that  do  not  indicate  transformation  or  amorphization  occurs 
in  boron  carbide  for  purely  hydrostatic  loading  [11,  12].  A  measure  of  intrinsic  stability  [6-8]  under 
hydrostatic  loading  requires  that  the  slope  of  the  pressure-volume  curve  remains  negative; 

-dp  dv  =  -dp  diA  V)  =  -V{dpldA){dAf  >  0  ^  dp/dA  <  0  (33) 

This  stability  condition,  which  correlates  with  a  positive  incremental/tangent  bulk  modulus,  is  clearly 
satisfied  by  the  solutions  in  Figure  1(a),  demonstrating  mechanical  stability  of  boron  carbide  under 
purely  hydrostatic  compression  to  large  pressures. 

3.2.  Simple  Shear  under  Constant  Compression 


Now  consider  spherical  deformation  of  an  element  of  initially  fully  crystalline  material  of  length 
L  superposed  with  shear  displacement  of  magnitude  u  that  potentially  varies  only  with  the  X 
coordinate.  In  Cartesian  coordinates,  the  deformation,  deformation  gradient  matrix,  Jacobian 
determinant,  and  shear  strain  are,  respectively, 

x  =  A^'^X  \a^'^  0  0  " 

y  =  A^'^Y  +  v(X),  F=  Y  0  ,  ^  =  detF>0,  Y  =  du/dX  (34) 

z  =  A^'^Z  0  0  A^^^ 

Here,  A  is  constant  with  respect  to  reference  coordinates  (VA  =  0),  but  y  may  generally  vary  with  X, 
where  Xe[0,  L].  Under  such  shearing  and  volumetric  deformation,  (18),  (19),  and  (21)  reduce  to 

F^(A,^)  =  F^,(A,^)  =  F^,(A,^) 

= HA/m] 

F^,{A,  Y,  ^)  =  {laoX(^)  -  koi(^)  ln[H/(^(^)]  (y  (35) 

P^,^(Y,^)  =  Kx(^)[(^(Or'^}Y  =  T 

PxZ  =  PyZ  =P.X=PzY=^ 


q(A,  Y,^)=<i>  'a)  {fio  X(0/c^(^)  -  HAlm'] } 

+  XX^H{f[c^(^)r''^'''-f-ln[^/(^(^)]  +  iY^}  +  iiX^)^o{ln[^/(^(^^ 


(36) 


p(A,j,^)  =  M^)/A]  {1  -  [Almr  -  \ }  -  [^oi(0/^]  ln[^/(^(^)]  (37) 

In  the  third  line  of  (35),  shear  stress  %  is  work  conjugate  to  y  in  the  energy  increment  per  unit 
reference  volume  P\dF.  Stress  equilibrium  (8)  with  (34)  and  (35),  presuming  stresses,  like 
deformation,  can  only  vary  with  the  X  coordinate,  reduces  to 

dF^^  ldX  =  dF^^ldX  =  d%ldX  =  d  (38) 

Order  parameter  equilibrium  in  the  first  of  (9)  with  (15)  and  (36)  becomes 

5=-U'/r)f[n„x/<l>-HoXr”"-4'"  -X„(i/(^)ln(.4/(|))] 

The  same  interpolation  functions  first  introduced  in  (30)-(32)  are  invoked  the  remainder  of  this  work. 
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First  sought  are  spatially  homogeneous  solutions  to  (35)-(39)  wherein  deformation  gradient  and 
order  parameter  do  not  vary  with  position  X: 

A  =  constant,  y  =  constant,  ^  =  constant  (40) 


From  (35),  such  solutions  clearly  satisfy  stress  equilibrium  (38).  Also  compared  later  is  the  ratio  of 
total  energy  with  possible  transformation  to  elastic  energy  for  the  case  when  no  phase  transformation 
occurs,  where  for  the  spatially  homogeneous  solution. 


A>  =  A>/V  =  (1/F)  //]dF  =  W{A,^,^)  +  Ti,^ll 

W,(A,y)  =  i^,[3(A^^^-l)-2lnA  +  y^]  +  jl,(lnAf 


Figure  2.  Shear  stress  (a),  order  parameter  (b),  Cauchy  pressure  (c),  and  ratio  of 
total  energy  to  energy  for  purely  elastic  deformation  (d)  for  simple  shear  y  under 
spherical  compression  A. 

Homogeneous  solution  data  are  obtained  as  follows  for  shear  and  compressive  loading.  Volume 
is  reduced  incrementally  by  reducing  A  from  unity  in  decrements  of  0.02.  For  each  volume 
decrement,  shear  y  is  increased  incrementally  from  zero  to  0.15.  Equilibrium  equation  (39)  is  solved 
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via  simple  numerical  iteration  to  obtain  order  parameter  Once  ^  is  known,  functions  (30)-(32)  are 
evaluated  and  then  substituted  into  analytical  expressions  (35),  (37),  and  (41)  to  determine  stress 
components,  pressure,  and  energies.  Shown  in  Figure  2(a)  is  the  shear  stress  t  (normalized  by  initial 
shear  modulus)  predicted  from  the  phase  field  model  at  various  compressive  strains  A. 
Complementary  results  for  the  order  parameter  and  pressure  (normalized  by  initial  bulk  modulus)  are 
given  in  Figures  2(b)  and  2(c).  As  volumetric  compression  increases  (i.e.,  as  A  is  reduced),  pressure 
increases,  transformation  to  the  amorphous  phase  becomes  more  rapid,  and  shear  strength  diminishes. 
Complete  amorphization  is  approached  as  shear  strains  become  large,  i.e.,  ^^1  for  y>0.15.  Such 
phase  field  predictions  agree  with  trends  observed  in  DFT  [11,  12]  that  amorphization  in  boron 
carbide  is  triggered  by  non-hydrostatic  loading  involving  shear  strain  and  is  accelerated  by 
increasing  superposed  compressive  stress.  As  indicated  by  ratio  values  less  than  unity  in  Figure  2(d), 
the  total  energy  of  the  phase  field  solution  is  always  less  than  that  of  the  purely  elastic  solution  (^  =  0) 
for  neo-Hookean  elasticity,  demonstrating  metastability  of  the  latter  model  relative  to  the  phase  field 
model.  Intrinsic  stability  at  constant  volumetric  deformation  A  under  dead  loading  by  shear  stress  t 
requires  [6,  8,  16] 

dP :  dF  =  dP^^dF,,;^,  =  dr  dy  =  (0T/fiy)(dy)^  >  0  =>  dx/dy  >  0  (42) 

This  condition,  which  correlates  with  a  positive  incremental/tangent  shear  modulus,  is  first  violated 
by  each  of  the  solutions  in  Figure  2(a)  when  shear  stress  t  reaches  a  maximum  or  critical  value  tc, 
demonstrating  mechanical  instability  of  boron  carbide  under  simple  shear  loading  possibly 
superposed  on  constant  volumetric  compression.  Peak  stresses  at  the  onset  of  such  instability  are 
listed  along  with  other  solution  variables  in  Table  2.  As  volume  decreases,  peak  shear  stress 
decreases,  pressure  increases,  and  the  value  of  the  order  parameter  at  instability  increases. 


Table  2.  Homogeneous  phase  Held  solution  at  critical  loading  or  instability  point. 


Volume  ratio  A 

Peak  stress  tc  [GPa] 

Shear  yc 

Order  parameter 

Pressure  pc  [GPa] 

1.00 

3.91 

0.36 

-3.02 

0.98 

2.11 

0.53 

-0.43 

0.96 

0.84 

0.046 

0.70 

2.33 

0.94 

0.28 

0.052 

0.84 

6.35 

0.92 

0.12 

0.067 

0.91 

11.78 

0.90 

0.07 

0.081 

0.94 

17.96 

Now  sought  are  spatially  heterogeneous  =  ^(X)]  but  stress-free  solutions  corresponding  to 
vanishing  of  P  in  (35)  under  null  compression  and  omitting  volume  collapse  upon  transformation: 

A  =  E^l^  det  F  =  (j)  =  1 
^PxX=  PyY  =Pzz=^ 

<1,  =  P.Y  =  Pyx  (Y.  ^)  =  PoX(Oy  =  Po  (1  -  Y  =  0 

q(y,^)  =  iPoX’(OY^=-Poa-^)Y^ 

The  final  two  lines  in  (43)  combine  to  imply  elastic  driving  force  in  (36)  must  also  vanish  and  that 
shear  strain  can  be  nonzero  only  in  regions  where  amorphization  is  complete: 
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=  0,  5| 


X=L12 


=  1 


^W  =  - 


(45) 


(46) 


q  =  -TY/(l-^)  =  0,  Y(X)  =  0V^(X)<1  (44) 

Order  parameter  equilibrium  in  the  first  of  (9)  becomes,  for  the  present  one-dimensional  shear 
problem  with  r  =  0  boundary  conditions  and  full  amorphization  at  the  midpoint  of  the  domain, 

dXjdX^  -  ^le  =  0 ,  {d^ldX)l__^  =  {d^ldX)\ 

This  is  a  homogeneous,  ordinary  second-order  differential  equation  with  exact  solution 
[ [£^42/)  ^  g-i/(2/)  yi  ^  ^-xn ^  for  X  <  Z/2 

^  ^3./(2/)  ^-1  ^xn  ^  ^  ^-1  ^-xn  x>Lll 

Total  energy  per  unit  cross  sectional  area  is,  noting  that  (43)  and  (44)  imply  W=  0, 

f  r[^7^  +  iWdXf  ]dX  =  2r(e^''  -  »  2r  (47) 

where  the  final  approximation  effectively  becomes  exact  for  L>  10/.  This  calculation  verifies  that  T 
is  indeed  the  surface  energy  of  a  localized  amorphous  band,  the  factor  of  two  in  (47)  accounting  for 
each  side  of  the  band  surface.  Shown  in  Figure  3(a)  is  the  heterogeneous  solution  (46)  for  a  domain 
length  of  T  =  25/,  verifying  that  the  width  of  the  localized  amorphous  band  wherein  ^  >  0.75  is 
approximately  /.  Shown  in  Figure  3(b)  is  the  ratio  of  total  energy  from  the  homogeneous  phase  field 
solution  under  null  volumetric  compression  (i.e.,  pure  simple  shear  loading)  to  the  total  surface 
energy  of  the  localized  solution  verified  in  (47).  When  this  ratio  in  Figure  3(b)  exceeds  a  value  of 
unity — ^which  occurs  at  y  »  0.0125 — the  localized  stress-free  solution  is  energetically  favorable  to  the 
homogeneous  solution,  suggesting  local  transformation  and  fracture  may  be  expected. 


Figure  3.  Order  parameter  profile  for  localized  stress-free  solution  (a)  and  ratio  of 
total  energy  for  homogeneous  phase  field  solution  to  surface  energy  (h)  under 
simple  shear  y- 


4.  Discussion 


The  phase  field  theory  developed  and  applied  in  the  present  paper  consists  of  the  following 
principal  components; 

•  Order  parameter  denoting  transformation  from  crystal  to  amorphous/glass  phase  (2) 
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•  Volumetric  inelastic  deformation  (11)  associated  with  mass  density  increase  upon  amorphization 

•  Compressible  neo-Hookean  nonlinear  elastic  strain  energy  (13) 

•  Interpolation  function  (30)  relating  order  parameter  increase  and  volume  decrease  (spherical 
eigenstrain)  similar  to  that  used  to  interpolate  eigenshear  in  twinning  models  [19-21] 

•  Degraded  shear  modulus  and  fixed  bulk  modulus  upon  amorphization 

•  Interpolation  function  (31)  relating  order  parameter  increase  and  shear  modulus  decrease  similar 
to  that  used  to  reduce  elastic  modulus  in  fracture  models  [22-25,  30] 

The  only  parameters  entering  the  phase  field  component  of  the  theory  are  surface  energy  and 
localization  width,  which  are  obtained  directly  from  atomic  theory  [9]  and  experiment  [4,  14],  The 
initial  bulk  modulus  of  the  neo-Hookean  potential  is  fit  to  uniquely  to  DFT  pressure-volume  data 
[11],  and  the  initial  shear  modulus  is  obtained  directly  from  second-order  elastic  constants  of  DFT 
[10,  11].  Thus,  the  phase  field  model  is  considered  purely  predictive,  with  no  adjustable  parameters. 
Of  course,  interpolation  functions  such  as  (30)-(32)  could  be  generalized  or  altered  as  needed  in  the 
future  to  match  additional  data  for  other  stress-deformation  states  as  such  data  become  available 
from  more  extensive  experimentation  and  atomic  simulation.  Furthermore,  the  model  should  be 
implemented  in  a  numerical  framework  such  as  the  finite  element  method  [19-21,  25]  to  study 
boundary  value  problems  involving  more  complicated  geometries  and  boundary  conditions. 
Extension  of  the  theory  towards  a  homogenized  polycrystal  model  [16,  44]  wherein  multiple 
anisotropic  B4C  crystals  occupy  a  single  material  point  in  a  computation  is  also  foreseeable. 

The  present  results  suggest  that  boron  carbide,  under  simple  shear  loading,  should  become 
unstable,  undergo  substantial  amorphization,  and  demonstrate  severe  strength  loss  at  shear  strains  on 
the  order  of  5%,  with  peak  shear  strengths  on  the  order  of  4  GPa,  the  latter  value  decreasing  with 
increasing  superposed  pressure  (Table  2).  While  the  trends  in  phase  field  results  agree  with  DFT 
predictions,  maximum  shear  stresses  observed  in  DFT  can  exceed  20  GPa,  with  shear  strains  at 
instability  or  peak  load  on  the  order  of  10-15  %  [11].  It  is  suggested  that  the  present  phase  field 
predictions  are  more  physically  realistic  when  compared  to  real  solids  with  defects,  since  it  is 
unlikely  that  any  brittle  material  could  sustain  such  large  shear  stresses  without  undergoing  mode  II 
fracture  which  appears  prohibited  by  the  small  system  size  and  boundary  conditions  used  in  DFT.  A 
complementary  interpretation  is  that  the  present  phase  field  model  of  shear  degradation  accounts  for 
simultaneous  amorphization  and  fracture  as  suggested  by  the  form  of  (31),  whereas  prior  results 
reported  from  DFT  [4,  10-12]  account  for  structure  collapse  but  not  complete  bond  breaking. 

The  present  theory  and  solutions  do  not  account  explicitly  for  the  presence  of  defects  in  the 
crystal,  such  as  point  defects,  dislocations,  stacking  faults,  and  twin  boundaries,  all  of  which  are 
known  to  exist  in  boron  carbide  ceramics  [4,  14,  32,  33].  The  model  and  analysis  may  be  altered  in  a 
simple  empirical  way  to  account  for  defects  simply  by  increasing  or  decreasing  the  surface  energy  F 
according  to  whether  initial  defect  concentrations  increase  or  decrease  the  energetic  barrier  for 
localization.  For  example,  if  a  stacking  fault  promotes  local  amorphization,  F  may  be  lowered  by 
some  suitable  fraction  of  the  stacking  fault  energy.  Alternatively,  for  a  more  realistic  and  predictive 
representation  of  effects  of  defects,  the  phase  field  theory  should  be  extended  to  include  additional 
order  parameters  describing  (partial)  dislocations  and  deformation  twins  [18-21,  45],  and  numerical 
simulations  wherein  initial  defect  geometries  are  resolved  explicitly  should  be  considered,  as  has 
been  done  elsewhere  in  a  study  of  heterogeneous  stress  distributions  in  highly  twinned,  slip-cast 
boron  carbide  [33].  Regardless,  dislocation-type  defects  are  thought  to  be  relatively  immobile  as 
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evidenced  by  the  low  ductility  and  high  stiffness  of  boron  carbide. 

The  model  and  analysis  reported  in  this  work  has  been  limited  to  quasi-static  loading.  In  impact 
or  ballistic  loading,  local  stress  field  perturbations  associated  with  elastic  wave  interactions  (e.g., 
shear  and  bulk  compression  waves)  and  heterogeneities  (e.g.,  crystal  defects  and  grain  boundaries) 
are  expected  to  promote  amorphization.  This  expectation  is  supported  by  previous  nonlinear  analysis 
and  simulations  invoking  a  Born-like  instability  criterion  for  localization  [16,  17],  wherein  dynamic 
deformation  facilitated  instability  relative  to  static  loading.  The  level  of  amorphization  or  structural 
disorder  in  boron  carbide  has  also  been  reported  as  greater  in  dynamic  indentation  experiments 
compared  to  static  indentation  experiments  [13,  46]. 

5.  Conclusion 

A  nonlinear  phase  field  theory  has  been  developed  and  applied  towards  amorphization  and  shear 
failure  in  boron  carbide  ceramic.  Analytical  and  iterative  numerical  solutions  have  been  obtained  for 
hydrostatic  loading  and  for  simple  shear  loading  in  conjunction  with  constant  volumetric 
compression.  Model  predictions  agree  with  trends  reported  elsewhere  in  experiments  and  quantum 
mechanical  (DFT)  simulations:  transformation  does  not  occur  under  purely  hydrostatic  loading; 
transformation  occurs  under  shear  loading  and  is  promoted  by  superposed  compressive  stress;  the 
tangent  shear  modulus  and  shear  strength  decrease  and  the  material  response  becomes  unstable  upon 
application  of  sufficient  shear  stress  or  shear  strain.  The  proposed  theory,  which  essentially  includes 
no  adjustable  or  free  parameters,  appears  to  be  the  first  phase  field  model  of  pressure-shear  induced 
amorphization  in  crystalline  solids. 
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